Charmonium spectral functions with the variational method 
in zero and finite temperature lattice QCD 
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We propose a method to evaluate spectral functions on the lattice based on a variational method. 
On a lattice with a finite spatial extent, spectral functions consist of discrete spectra only. Adopting 
a variational method, we calculate the locations and the heights of spectral functions at low-lying 
discrete spectra. We first test the method in the case of analytically solvable free Wilson quarks 
at zero and finite temperatures and confirm that the method well reproduces the analytic results 
for low-lying spectra. We find that we can systematically improve the results by increasing the 
number of trial states. We then apply the method to calculate the charmonium spectral functions 
for S and P-wave states at zero-temperature in quenched QCD and compare the results with those 
obtained using the conventional maximum entropy method (MEM). The results for the ground state 
are consistent with the location and the area of the first peak in spectral functions from the MEM, 
while the variational method leads to a mass which is closer to the experimental value for the first 
excited state. We also investigate the temperature dependence of the spectral functions for S-wave 
states below and above Tc- We obtain no clear evidences for dissociation of J/i/" and rjc up to lATc- 

PACS numbers: 11.15.Ha,12.38.Gc,12.38.Mh 



I. INTRODUCTION 

Quark-gluon-plasina (QGP) is expected to be formed 
at sufficiently high temperatures and densities and con- 
sidered to play an important role in the early universe 
and in the core of neutron stars. In heavy-ion collision 
experiments, the QGP created is expected to cause a 
suppression of J/ip creation due to the dissociation of 
charmonia in the deconfinement phase of QCD Q. Sup- 
pression of J/ ip particles was actually observed in several 
experiments at the Super Proton Synchrotron and the 
Relativistic Heavy Ion Collider [s^. Here, not only J/iJj 
but also Xc and ip' states contribute to the total yield of 
J/tp Q. To investigate properties of these charmonia in 
the deconfinement phase, hadronic spectral functions are 
studied at zero and finite temperatures. 

spectral functions on the lattice are calculated conven- 
tionally by using the maximum entropy method (MEM) 
in which spectral functions are extracted from Eu- 
clidean correlation functions using a Bayesian probabil- 
ity theory. With the MEM, the temperature dependence 
of charmonia spectral functions has been investigated in 
quenched [E-Q and two- flavor QCD From these 

studies, the S-wave charmonia {rjc, J/ip) are suggested 
to survive up to temperatures higher than about l.STc, 
where Tc is the critical temperature. On the other hand, 
the P-wave charmonia {xco, Xci) are suggested to dissolve 
just above 

The spectral functions from the MEM are continuous 
as expected in infinite- volume field theories. In the finite 
volume, on the other hand, the spectral functions consist 
of discrete spectra only, since the degrees of freedom of 



the theory is finite. Therefore, the meaning of the con- 
tinuous spectral functions from the MEM is not quite 
clear. This is problematic, in particular, in cases when 
the MEM leads to ambiguous results depending on, e.g., 
the choice of the default model. 

In this paper, we propose a method to directly calcu- 
late spectral functions at discrete spectra on the lattice, 
applying a variational method. Variational method is a 
powerful tool to extract information of low-lying discrete 
spectra (T]| . In our previous paper, we studied the tem- 
perature dependence of charmonium wave functions for 
the ground and first excited states using a variational 
method [l^]. We extend the method to evaluate spec- 
tral functions, i.e., the location and the height for each 
low-lying discrete spectra. 

With the variational method, information for low-lying 
spectra is extracted by diagonalizing correlation matri- 
ces between various smeared operators. The spectral 
functions are defined in terms of correlation functions 
between pointlike operators. We thus include pointlike 
operators among the smeared operators. We calculate 
spectral functions from the element of the correlation 
matrices corresponding to the pointlike source and sink 
operators. We show that the accuracy and reliability of 
the resulting spectral functions can be systematically im- 
proved by increasing the number of smeared operators. 

At nonzero temperatures, the spectral function for a 
meson operator Or is given by 
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where Ek is the energy of the state |fc) and Z{T) is the 
finite temperature partition function. As T is increased, 
we expect that the height of each pole varies and, simul- 
taneously, more poles at Ek — Ek' appear in addition to 
the zero-temperature poles at Ek — Eq, where Eq is the 
ground-state energy. When a particle dissociates, we may 
expect that the corresponding spectral function shows a 
drastic change around the dissociation temperature: The 
height of the pole corresponding to that particle will be- 
come lower and will be eventually buried by nearby poles. 
Our final goal is to clarify the fate of various charmonia 
at high temperatures. 

This paper is organized as follows. We introduce our 
method to calculate locations and heights of the peaks 
of meson spectral functions with the variational method 
in the next section. In Sec. IIIIl we first test the method 
for the case of free Wilson quarks. We construct me- 
son correlation matrices using various Gaussian smearing 
functions and compare the results of spectral functions 
with analytic solutions. We then apply the method to 
calculate the charmonium spectral functions in quenched 
QCD at zero and finite temperatures in Sec. IIVI Our 
conclusions are given in Sec. |Vl A preliminary report 
was presented in [l^. 



II. SPECTRAL FUNCTIONS WITH THE 
VARIATIONAL METHOD 

In this section, we introduce a method to calculate lo- 
cations and heights of the peaks of meson spectral func- 
tions with a variational method [TT'| . 

The Euclidean meson correlation function is defined by 



CT{x,t) = (Or(2r,t)OM0,0)), 



(2) 



where OY{x,t) = q{x,t)Tq{x,t) is the pointlike meson 
operator for channel F and F = 75, 7^, 1, 757^ (j = 1, 2, 3) 
correspond to pseudoscalar (Ps), vector (Ve), scalar (Sc) 
and axial-vector (Av) channels, respectively. For Ve and 
Av channels, we average the correlation functions over 
i = 1,2,3. Then, the meson spectral function pr{uj,p) is 
related to its Fourier transform 



Cr{p,t) = / d^xCr{x,t) e'P''- 



(3) 



by 



^ f°° , ^ cosh[cj(i - ^)] 
Jo smh[^\ 

where T is the temperature. Setting the lattice spacing 
a = 1, we have 1/T = Nt with Nt the temporal lattice 
size. 

On finite lattices, the spectral function consists of dis- 
crete spectra only. Therefore, the integral over w in 
is actually a summation over discrete values of u. For 



the zero momentum case, p = 0, we rewrite (j4]) on finite 
lattices as 

r, sr ( \ cosh[mfcft- A^t/2)] 

C^r W = V PrKruk) — —— — 5 

^ smh[TOfeiVt/2J 

with fc = 1, 2, • • •, where /5r(w, 0) = pT{mk)5{ui — mk) ■ 



TABLE I: Smearing parameters Ai for the test of the varia- 
tional method. A\ corresponds to the point operator. 
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FIG. 1: t- and n-dependence of m| (top) and pT(m1 ) (bot- 
tom) obtained on the 20^ x 128 lattice for the Ps channel. 



Let us introduce smeared meson operators 

Ot{x, t)i = ^ LUi{y) uj.i{z) q{x + y,t)r q{x + z, t) (6) 



with an appropriate gauge fixing j22| . where uji{x) {i ~ 
1, 2, • • • , n) are smearing functions. We define an n x n 
meson correlation matrix Cr(t) = [C'r(i)i.j] by 



Cr(% =^(Or(f,0.(^r(0,0), 



(7) 



In this study, we set wi(x) — S{x), so that Cr(i)ii is only 
the Cr{t) defined by (0). 

By solving a generalized eigenvalue problem 



Cr(t)vW = Afc(t;to)Cr(to)v 



(fc) 



(8) 
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FIG. 2: The same as Fig. [T] for the Sc channel. 



for fc = 1, • 

Afc(t; to) 



n, we define efi'ective masses m'jf (t; to) by 
cosh[mf{t;to){t-Nt/2)] 



cosh[mf{t;to){to - Nt/2)y 



(9) 



Denoting A = diag{Ai, • • • , A„} and V = [v^^) • • • v(")], 
we rewrite dH) as Cr{t) = Cr(to)VAV-i. Then, the 
(1,1) element of this relation reads 



CrWii 

^ cosh[mf{t;to){to ~ Nt/2)] 

cosh[mfit;to)it-Nt/2)] 



smh[mfit;to)Nt/2] 



8mh[mf{t;to)Nt/2] 



(10) 



Comparing (O and pH)). we define an effective spectral 
function 
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FIG. 3: Meson spectral function at three lowest-lying states 
in the Ps and Ve channels obtained at t = 63 on a 20^ X 128 
lattice in the free quark case. Cross, square, circle, triangle, 
and downward triangle symbols are the results of the varia- 
tional method for n = 3,4, 5, 6, 7, respectively. Symbols in the 
bracket mean that the corresponding signals are not asymp- 
totic even at t = 63. The analytic solutions are shown by the 
asterisks. 



In the calculations shown below, we further apply the 
midpoint subtraction method [l^ to the meson corre- 
lator matrices in order to subtract the constant mode 
contributions: 



pr{mf{tM) = (Cr(io)V),,(V-i)fei 

^ sinh[mf{t;to)Nt/2] 

cosh[mf{t;tQ){to-Nt/2)y^ ' 

When we let t and sufficiently large, toS.^ approaches 
to the mass of the fc-th state nik and pri'^k) approaches 
to the spectral function priiTik) defined by ([5]). 

Note that while only the (1,1) element is related to 
pr(TOfc*), all n trial states contribute in (fTT|) . Keeping n 
finite introduces a systematic error in the location and 
height of low-lying spectra at finite t. By increasing n, 
we can systematically improve the results. On the other 
hand, we note that setting n too large can lead to large 
statistical fluctuations and numerical instabilities. 



Cr{t) ^ Cr{t)-Cr{Nt/2). (12) 

Accordingly, cosh in ([5]) and (|9))- (ITT]) should be modified 
as cosh [m{t - Nt/2)] -> cosh [m{t ~ Nt/2)] - 1 while 
other factors including the sinh terms remain unchanged. 
In principle, we can treat the constant mode as an eigen- 
state in the variational method too. We have confirmed 
that the low-lying physical modes from variational calcu- 
lations with and without the midpoint subtraction pro- 
cedure are consistent with each other. We find, however, 
that the midpoint subtraction improves the arithmetic 
precision of signals at large t and thus the resulting effec- 
tive spectral functions are more stable with the midpoint 
subtraction. 
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FIG. 4: The same as FIG. [3] for the Sc and Av channels. 



III. TEST WITH FREE QUARKS 

In order to test the method, we first study the case 
of free Wilson quarks. We compare spectral functions 
from the variational method with the analytic solutions 
on anisotropic lattices with the anisotropy ^ = as /at = 4, 
where and at are the spatial and temporal lattice spac- 
ings, respectively. At the Wilson parameter r = 1, we 
adjust the quark mass to « 0.7501 to approximately re- 
produce the grand-state meson masses in quenched QCD 
studied in the next section. 

The analytic solutions for meson spectral functions 
with free Wilson quarks are given in Appendix |^ for 
the Ps, Ve, Sc and Av channels. The spectral functions 
with the variational method are calculated by setting link 
variables to unity in the QCD code to be used in the next 
section. In this paper, we adopt Gaussian smearing func- 
tions defined by 
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FIG. 5: t- and n-dependence of rn^ (top) and pr(m^^) (bot- 
tom) obtained on the 20^ x 32 lattice for the Ps channel. 
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FIG. 6: The same as Fig. [5] for the Sc channel. 



z = l,2,. 



(13) 



with the smearing parameters Ai listed in Table U] Here 
Ai = oo is for the point operator. 

We first test on a 20^ x 128 lattice with ^ = 4. In Figs.[T] 
and[21 we show the effective mass m'^{t, to) (top) and the 
effective spectral function pr{m'i^{t, to)) (bottom) for the 
lowest three states in the Ps and Sc channels as functions 



of t. We extract the signals at the largest i = 63 (f = 64 
is excluded by the midpoint subtraction procedure). In 
order to suppress the contamination of higher states, we 
choose large to (< under the condition that the signals 
of and p{mf^) for the lowest three states are stable 
for all values of f up to n = 7 (see Appendix IB|) . The 
values of tg are given in Appendix [BJ The results for 
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FIG. 7: Meson spectral function at three lowest-lying states 
in the Ps channel obtained by the variational method with 
n = 7 at t = 15 on the 20^ x 32 lattice (left), and at f = 63 
on the 20"^ x 128 lattice (right) in the free quark case. The 
asterisk symbol is for the analytic solutions. 
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FIG. 8: The same as Fig. [7] for the Sc channel. 



the Ve (Av) channel are similar to those of the Ps (Sc) 
channel. 

From these figures, we find that the ground-state sig- 
nals mi and p{mi) can be safely extracted even with a 
small n. A larger n is required for excited states to ob- 
tain asymptotic signals at t = 63, in particular, for the 
second excited state (fc ~ 3) of P-waves. On the other 
hand, we find that setting n too large cause numerical 
instabilities due to the limitation of the arithmetic pre- 
cision - in the present test, n > 7 induces instabilities at 
large t for several states. We thus restrict ourselves to 
n<7. 

Figures |3] and |4] show the results of and pr(™fc) for 
the lowest three states obtained at t = 63 with n — 3, 
4, • • •, 7. The results in the bracket are those appar- 
ently not asymptotic at < = 63 from the t-dependence 
of the effective mass or the effective spectral function. 
The analytic solutions are given by the asterisks. We 
find that, except for the case of the second excited states 
in P-wave, the analytic solutions are well reproduced by 



TABLE II: The simulation parameters for the plaquette gauge 
action and the C'(a)-improved Wilson quark action on our 
anisotropic lattice. 
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choosing a sufficiently large n. For the second excited 
states in P-wave, we observe slight deviations from the 
analytic results even with our largest n (see Fig. S]). The 
fact that the deviations become smaller with increasing 
n suggests that these results are not fully asymptotic 
yet. A possible cause will be our choice of the trial states 
based on the Gaussian smearing functions which may not 
be overlapping well with P-wave excited states. Another 
reason may be the large contamination of the constant 
mode in P-wave correlation functions (see Appendix . 
Although the constant mode is removed by the midpoint 
subtraction procedure, the resulting signal suffers from 
lower precision. We leave these issues for future investi- 
gations. 

On finite temperature lattices with a small temporal 
extent Nt, the range of t available for the variational 
analyses is limited. To study its influences, we repeat 
the test on an anisotropic 20'^ x 32 lattice with ^ = 4, 
adopting the same simulation parameters. In the free 
quark case, because no interactions with the thermal 
background medium exist, we have no additional poles 
at T > in dH). 

Effective masses and effective spectral functions ob- 
tained on the finite temperature lattice are shown in 
Figs. [5] and IH] for the Ps and Sc channels. We adopt 
to = 14 and t = 15. In Figs. [7] and [51 we compare spectral 
functions on 20^^ x 32 and 20'^ x 128 lattices for the lowest 
three states in these channels obtained with n = 7. The 
asterisks represent the analytic solutions. Results for the 
Ve and Av channels are similar to those of the Ps and 
Sc channels, respectively. We find that the limitation 
of the range of t requires a larger n to extract asymp- 
totic signals. On the Nt = 32 lattice, the second excited 
states show deviations from the analytic results even with 
n = 7. To overcome the problem, a set of more optimally 
smeared operators will be needed. By examining both t 
and n dependences of the results, however, we can get an 
idea to which extent the results are asymptotic. 



IV. CHARMONIUM SPECTRAL FUNCTIONS 

In this section, we study charmonium spectral func- 
tions in quenched QCD and compare the results of the 
variational method with those obtained by the conven- 
tional MEM. 
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FIG. 9: Charmonium spectral function calculated by MEM 
for the Ps and Ve channels at zero temperature. We use p4p 
as the default model function. The results for a — 0.5, 1.0, 2.0 
are shown by solid, dashed, and dotted lines, respectively. 



A. Simulation parameters 

We perform simulations on anisotropic 20^ x Nt lat- 
tices with the renormalized anisotropy ^ = 4 adopt- 
ing the standard plaquette gauge action. We study at 
(3 — 6.10 where the spatial lattic e sp acing determined 
by the Sommer scale tq = 0.5 fm [ill is = 0.0970(5) 
fm (a^^ = 2.030(13) GeV). Our spatial volume is thus 
about (2 fm)'^. The bare anisotropy for ^ = 4 is 
7G = 3.2108 |17j. For valence quarks, we adopt an 0(a)- 
improved Wilson quark action. We set the bare fermionic 
anisotropy "fp = 4.94 and tree-level tadpole improved 
clover coefficients ce — 3.164 and cb — 1.911 to realize 
C = 4 (see Ref. [H for the definitions of the coupling 
parameters). In this study, we set the Wilson parameter 
r = 1 to suppress lattice artifacts in excited charmonia 
[T^ . We study at k = 0.10109, which corresponds to the 
physical charm quark mass on an isotropic lattice with a 
similar spatial lattice spacing. 

For the temporal lattice size, we adopt Nt = 160 for the 
zero-temperature simulation, and Nt = 32, 26, and 20 for 
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FIG. 10: The same as Fig. [9] for the So and Av channels. 



finite-temperature simulations at T « 0.88Tc, l.lTc, and 
1.4Tc, respectively. Here the critical temperature Tc is 
determined by the peak position of Polyakov loop suscep- 
tibility and corresponds to Nt ^ 28. After 20 000 sweeps 
for thermalization, we generate 299 configurations at zero 
temperature and 800 configurations at finite tempera- 
tures separated by 500 pseudo-heat-bath sweeps. Our 
simulation parameters are summarized in Table |lll To 
calculate smeared operators, we use the Coulomb gauge. 
Statistical errors for spectral functions are estimated by 
a jackknife method. 



B. Zero temperature 

At zero temperature, we calculate the locations and 
the heights of the peaks for charmonia spectral functions 
up to the first excited state for Ps, Ve, Sc and Av chan- 
nels. We first calculate the spectral functions with the 
conventional MEM Q using meson correlation functions 
for point operators. We use the range t = 1-60 for Ps 
and Ve channels, and t = 3-60 for Sc and Av channels, 
because the latter correlation functions suffer from lat- 
tice artifacts at t '-^ 1. For the default model m{uj), we 
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FIG. 11: mf(t,to) (left-top), mf{t,to) (right-top), 
Pr{mf{t,to)) (left-bottom) and pr{mf{t,to)) (right- 
bottom) of chamonium for Ps channel at zero temperature 
with the variational method. The reference point is chosen at 
to — 5. Fit ranges to determine the locations and heights of 
the spectral peaks for charmonia are also shown by horizontal 
solid lines. 



adopt 



m{ijj) — arriDyi uP' 



(14) 



where muM = 4.2 for Ps and Sc channels and 2.4 for Ve 
and Av channels, as determined by the asymptotic be- 
havior of meson correlation functions in the perturbation 
theory at a = 1 [1, Q . We estimate statistical errors by 
the jackknife method. To estimate systematic errors due 
to the choice of the default model, we multiply the factor 
a and vary it in the range 0.5-2.0. We have checked that 
when the default model is fixed, the results are stable 
under variations of various parameters in the MEM. 

Our results of spectral functions by MEM for Ps, Ve, 
Sc and Av channels are shown in Figs. [5] and [TU] for the 
cases a = 1, 0.5 and 2.0. We find that for the S- waves, the 
peaks are well isolated up to the first excited states and 
the results are approximately stable. But, for P-waves, 
the peaks are not isolated and peaks for excited states are 
not stable. We identify with the peak position defined 
by the maxima of the spectral function, and primk) with 
the area of the peak at nik- For the P-waves for which 
the peaks are not well isolated, we divide the spectral 
functions into each peak at the minima to compute the 
area. 

To calculate spectral functions with the variational 
method, we adopt the Gaussian smearing function (jl3p 



FIG. 12: The same as FIG. [TT] for the Ve channel. 



with the smearing parameters listed in Table HI The ef- 
fective masses and effective spectral functions with Iq = 5 
are shown in Figs. [TTV[T4l To extract the plateau, we fit 
the data in the range shown by horizontal solid lines in 
Figs. [Tn [T4l for each charmonium. Here, we choose the 
same fit ranges for effective masses and corresponding ef- 
fective spectral functions. The fit ranges [iminj^max] are 
chosen as follows. We first determine tmax as the up- 
per bound of the plateau of m'j.^ and pri'iTi'if ) for each 
k. Then, we study fmin-dependence of and pri'Ttk)- 
As expected from Figs. 11-14, when n is small, we some- 
times observe that nik keeps decreasing with increasing 
tmin up to imax. When n is sufficiently large, however, 
TOfe and pv{'mk) are stable in a wide range in the sense 
that iniin-dependences are smaller than, or at least com- 
parable with, the small statistical errors, though the sta- 
tistical errors become large when tmin becomes close to 
imax- Because a simple criterion to define a plateau re- 
gion is not available, in this study, we just choose tmin, 
where x^/dof for becomes closest to 1, and reject the 
data by adding brackets in the plots when a systematic 
tendency to deviate from the plateau is observed beyond 
statistical errors by increasing tmin towards tmax- With 
the present lattice size and statistical accuracy, we cannot 
exclude a mild decrease of rrik and pr^mk) by choosing 
the fit range at larger t. Therefore, our and pr("^fe) 
should be regarded as upper bounds for them. Resulting 
ranges of the value of x^/dof are 0.82-1.1 and 0.045- 
0.54 for mfe and pr{mk), respectively. We confirm that 
the variations of and pr{rnk) under a change of t^in 
to iniin ± 1 are less than 0.2% and 2%, respectively. 
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0.07 i 




FIG. 13: The same as FIG.[TT]for the Sc channel. 



FIG. 14: The same as FIG.[TT]for the Av channel. 



Our results for charmonium spectral functions are 
summarized in Figs. [TSUISl The errors for the varia- 
tional method's data are smaller than the symbols. The 
symbols in the brackets are not asymptotic. For the re- 
sults from MEM, statistical errors for nik and pr(wfc) are 
given for each values of a. We note that the statistical 
errors in and prinik) have strong positive correla- 
tion, i.e., the errors actually shape a thin oval inclined 
rightwards. The experimental values of corresponding 
charmonium masses with their errors are shown by verti- 
cal dashed lines fio'] . In this study, the charm quark mass 
is adjusted to approximately reproduce the experimental 
J/ijj mass. On the other hand, masses of r)c, etc., show 
slight deviations from experiment, in accordance with the 
previous observation with 0(a)-improved Wilson quarks 
that the charmonium hyperfine splitting is smaller than 
experiment in quenched QCD (20j . 

We find that the results for the ground states are well 
consistent with each other between the MEM and the 
variational method with all n studied. On the other 
hand, for the first excited states, the variational method 
leads to results discrepant from those of MEM. For the 
S-waves, we find that the discrepancy is beyond the er- 
rors estimated with the MEM and becomes larger with 
increasing n. We find that the results of the variational 
method converge to a point close to the experimental 
values. This suggests that the inclusion of higher states 
in the variational approach helps improve the signals for 
excited states. We also note that the results of MEM 
approximately corresponds to those of the variational 
method in the limit of small n, in accordance with the 



fact that the studies with MEM are based on the infor- 
mation of point-source correlation functions only. The 
results for P-wave excited states are similar, but the er- 
rors in the MEM results are quite large to draw a definite 
conclusion. 



C. Finite temperature 

Finally, we study the charmonia spectral functions at 
finite temperature. We calculate ruk and pr(™fc) with the 
variational method using the same smearing functions 
with n — 7 and to — 5. 

In Figs. and \^U[ we show the results of m^*^(t,to) 
and primf^ {t,to)) for the Ps and Ve channels at T = 0, 
0.88Tc, l.lTc, and lATc- We see no clear temperature 
dependence in ml^ (upper panels) up to 1.4Tc. On the 
other hand, pr{mf^) (lower panels) show a slight shift 
between the temperatures below and above Tc- How- 
ever, the temperature dependence is not a drastic one as 
expected when a particle dissociates. Thus, these results 
suggest that both rjc and J/ip survive up to lATc- Be- 
cause the effective masses and effective spectral functions 
are not asymptotic up to the largest t available, we do 
not attempt to fit a plateau. 

The results for the first excited state Ps channel are 
shown in Fig. [5T] We also found that results for the 
Ve channel look similar. We note that both the effec- 
tive mass and effective spectral function show strong T- 
dependence above Tc- The results for P- waves are also 
similar to those for the S-wave first excited states. Even 
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FIG. 15: Charmonium spectral function at the ground and 
the first excited states for Ps channel at zero temperature 
obtained on the 20'^ x 160 lattice. The reference point is at 
to = 5. Cross, square, circle, triangle, and downward triangle 
symbols indicate the data by the variational method with 
n = 3, 4, 5, 6, 7, respectively, and solid symbols indicate the 
MEM results. The vertical dashed lines indicate the range of 
experimental masses for the 77c(lS) and 77c(2S) mesons. 
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FIG. 16: The same as Fig. [T5]for the Ve channel. The vertical 
dashed lines indicate the range of experimental masses for the 
J/^p{lS) and i/'(2S) mesons. 



for the ground states, as shown in Fig.[22]for the Sc chan- 
nel, we find the strong T-dependence above Tc- This may 
suggest appearances of additional poles. However, since 
we cannot extract asymptotic signals from the present 
data [2^, it is difficult to draw a definite conclusion on 
the fate of corresponding charmonia above Tc. 
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FIG. 17; The same as Fig. [T5]for the Sc channel. The vertical 
dashed lines indicate the range of experimental mass for the 
Xco(lP) meson. 
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FIG. 18: The same as Fig.[T5]for the Av channel. The vertical 
dashed lines indicate the range of experimental mass for the 
Xci(lP) meson. 



V. CONCLUSIONS 

We introduced a method to calculate meson spectral 
functions with the variational method. We first con- 
firmed by a test in the free quark case that the method 
reproduces the analytic solutions well for several low- 
lying states when the number of trial operators, n, is 
sufficiently large. By introducing more trial operators, 
we can systematically improve the signal. On the other 
hand, a judicious choice of the trial operators is needed 
to obtain an asymptotic signal within the range of avail- 
able t, in particular, for highly excited states in P-waves. 
This imposes a severe limitation on the applicability of 
the method at high temperatures. A good feature of 
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respectively. 



the method is, however, that we can judge to which ex- 
tent the results are regarded as asymptotic by examining 
the t- and n-dependences of effective masses and effective 
spectral functions. 

We then adopted the variational method to calculate 
charmonium spectral functions in quenched QCD. Com- 
paring the results of the variational method with those 
of the conventional MEM at zero temperature, we found 
that the location and the area of the ground state peak 
by MEM are well reproduced by the variational method. 
For the first excited states, we find that the variational 
method leads to spectra much closer to the experimental 
ones. We note that the results of MEM approximately 
corresponds to those of the variational method in the 
limit of small n. 

We also studied the temperature dependence of spec- 
tral functions. We found that the effective masses for 
the S-wave ground states show no T-dependence up to 
lATc- Corresponding effective spectral functions show a 
slight shift between below and above Tc- The absence 
of a drastic T-dependence suggests that rjc and J/ip do 
not dissociate up to 1.4Tc. However, the limitation in the 
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FIG. 20: The same as Fig. [TS] for the Ve channel. 



range of t poses a severe constraint to extract asymptotic 
signals. 

In systems with finite volume, continuum spectra in 
the infinite-volume limit must break up into discrete 
spectra. In the case of free quarks, we have confirmed 
the appearance of expected discrete spectra in the spec- 
tral function for meson operators. In QCD at T = 0, we 
do not expect such additional discrete spectra around 
the ground state because they appear only above the 
two-particle threshold in the same channel. Accordingly, 
we did not observe them around the states we studied. 
When a meson dissociates at high temperature, we ex- 
pect that the corresponding discrete spectrum changes 
into a broad continuous peak in the spectral function 
at that temperature in the infinite-volume limit. On fi- 
nite lattices, we thus expect the appearance of additional 
discrete spectra around the original spectrum. We have 
observed strong T-dependence of the signal at intermedi- 
ate distances in several channels above Tc, which may be 
suggesting the appearance of these additional spectra in 
these channels. To draw a definite conclusion, however, 
we need asymptotic signals for these spectra. When such 
asymptotic signals become available, it is important to 
check the volume dependence and the overlap with the 
meson operators used. More work is needed, in particu- 



11 



0.56 
0.52 
0.48 

CNJ 

O) 

E 

0.44 
0.4 
0.36 

0.35 

0.3 

0.25 

jT^ 0.2 
E, 

k 0.15 
0.1 
0.05 




1 1 1 1 1 1 1 
- Ps, 1 ^' excited state 

Si IT 

■ * * 1 i 


1 1 1 1 1 
T=0 
0.88T. 
1.1T^ 


1 1 1 
1 — 1 — 1 

i U i 


* 1 

□ 


? 1 * * JK 

* ^ 3K x 


X _ 

* * + . 




[p 










1 1 1 



01 2345678 9 10 11 12 13 14 15 16 
t 



"T — I — I — I — I — I — I — I — I — I — I — I — I — r 



Ps, 1 excited state 



T=0 I — H 



0.88T. --X--- 
1.1T. 

1-4T, I B i 



□ 9 



f i I f 8 i I i f 



SK X k )i( ^ 



01 2345678 9 10 11 12 13 14 15 16 
t 



0.6 



0.5 



i 0.4 
E 



0.3 



0.2 



"T 1 1 1 1 1 1 1 1 1 1 1 1 1 T 

Sc, ground state 



T=0 I — I — I 
0.88T. --X— ■ 
1.1T. 

1-4T, ■: B i 



« « S 



□ 5 i * .. 



_i I I I I I I I I I I I I I ^ 



1 2 3 4 5 6 7 8 910111213141516 



0.1 
0.08 
0.06 

E, 

u 

0.04 
0.02 




"T — I — I — I — I — I — I — I — I — I — I — I — I — I — r 

7=0 I — I — I 
0.88T' --X— ' 



Sc, ground state 

f I I 
i ^ ^ I 



1-4^ I B ^ 



* ? I i u 



't'ttlfl I [t]jK)K)K) | <:)K 



01 2345678 9 10 11 12 13 14 15 16 
t 



FIG. 21: mf{t,to) and pr(mf{t,to)) for the Ps channel. FIG. 22: mf{t,to) and pT{mf{t,to)) for the Sc channel. 



lar, in optimizing the trial operators, to discuss the fate 
of charmonia above Tr. 
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Appendix A: Analytic solution of meson spectral 
functions for free Wilson quarks 

The free quark propagator for the Wilson quarks in the 
momentum space is given by 



S{P) 



—174 sinp4 — iV{p) + 1 — cosp4 + M-{p) 
sin2p4 + V^ip) + [1 - cosp4 + M{p)Y 



(Al) 



TABLE III: a, bj and c defined by for Ps, Ve, Sc, Av 

channels, where i is the spatial direction of F for Ve and Av 
channels. 





F 


a 


bj 


c 


Ps 


75 


1 








Ve 


li 


1 


kj 





Sc 


1 





1 


1 


Av 


75 7i 





1 - S^j 


1 



with 



1 

'Pip) = 



M{p) 



r ^^(1 — cospj) + m 



(A2) 



where ^ = Os/at is the lattice anisotropy _21]. With the 
antiperiodic boundary condition in the temporal direc- 
tion, the Fourier transform of the quark propagator in p4 
reads 
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3 = 1 

15m 



2(1 + X(p))' 



sin\y[E{p){t-Nt/2)] 
(A3) 



Wilson quark case is now given by 

~ ( \ 



N, ^ smh[E{p}Nt] 



where 

Su{P) 



(1 + cosh^[£;(p)iVt/2] 



sinh ^^(p) 



2(1 + sinh £;(p) cosh[£;(p)iVt/2] ' 

sinpj 

2(1 + M(p)) sinh £;(p5cosh[£;(p)iVf/2]' 
l-cosh£;(p) + A^(p) 



:^ sinh E{p) 



TABLE IV: to for the 20-^ x 128 lattice 



(A9) 



2{1 + M{p)) sinh E{p) cosh[E{p}Nt/2] 



Here, E{p) is the location of the pole of (jAl[) : 



cosh E{p) = 1 



2{1 + M{p)) 



(A4) 



(A5) 





n = 3 


4 


5 


6 


7 


Ps 


62 


62 


62 


60 


53 


Ve 


62 


62 


62 


56 


52 


Sc 


62 


62 


57 


25 


13 


Av 


62 


62 


57 


25 
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Thus, the meson correlation functions ([2]) in the free 
quark case are given by 



Crit) 



(A6) 



for flavor nonsinglet channels. Substituting (jASp . we ob- 
tain 



Crit) 



Nr. 



N! ? (1 + Mip)y cosh'^[E{p)Nt/2] 



1 



a — > bj — 



sin^ Pj 



3 



^cosh[2i?(pO(i-iVt/2)] 

smh E (p) 




sin^ Pj 



^^^'sinh2i?(p) 



(A7) 



for t > 0, where a, bj and c are calculated by traces of 7 
matrices as 

a = itr{rrt}-tr{r74rt74}], 
b, = i[tr{rrt} + tr{r7,rt7j]. 



;[tr{rrt} + tr{r74r^74}]. 



(A8) 



Concrete values of them are summarized in TABLE IIIII 
The last term in (jA7p is removed by the midpoint sub- 
traction procedure p^ . For the Ps channel, the last 
term is absent because 6j = c = as listed in Table Hill 
On the other hand, for the P-wave states (the Sc and Av 
channels), the last term is numerically large. 

Adopting the midpoint subtraction procedure, the an- 
alytic solution for the meson spectral function in the free 



FIG. 23: Effective mass for the Av channel on the 20^ x 128 
lattice with to = 40. All the seven levels for n = 7 are shown 
as functions of t. 



Appendix B: Effective meson masses and meson 
spectral functions in the free quark case 

As discussed in Sec. IIIII we choose to as large as possi- 
ble in the region where signals of and p{m'^) up to 
the second excited states (A; = 1, 2, 3) are both stable for 
all t and n up to n = 7. We use t — Nt/2 — 1 to extract 
asymptotic values of ruk and pr{mk). Our choices of to 
for the 20"^ x 128 lattice are summarized in Table ITVl On 
the 20'^ X 32 lattice, we adopt to = 14 for all channels 
and all values of n. 

Results of the effective mass mf^ and the effective spec- 
tral function pri'm')^) as functions of t and n are given 
in Figs. [T] and [5] for the Ps and Sc channels. Results for 
the Ve and Av channels are similar to those for the Ps 
and Sc channels, respectively. 



13 



When we adopt to larger than the value given in Ta- 
ble IIV[ we encounter unstable signals at several interme- 
diate values of t. An example is shown in Fig. [23l We 
find that a strange level appears at intermediate values 
of As we vary t, the strange level crosses the ordinary 
levels which have milder dependences on t. We find that 



the strange levels are suppressed when we adopt a suf- 
ficiently small to, or limit ourselves to smaller values of 
n. These strange levels disturb the naming of low-lying 
states and the reliability of the signals. To avoid such 
levels for all n up to n = 7, we adopt the values of to 
listed in Table HVl 
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